rm(list = ls())
library(foreign)
library(stargazer)
library(Amelia)

set.seed(981237)
source('functions.R')
source('pool_compare.R')
load('main_data.Rdata')

#----------------------------------------------------------------------------------------
#Footnote 1:
#----------------------------------------------------------------------------------------
mean(cus_spells[cus_spells$cus_sr==1, 'cus_t_surv'])
#Average duration of revolutionary regimes is 39.1 years
mean(cus_spells[cus_spells$cus_sr==0, 'cus_t_surv'])
#Average duration of non-revolutionary regimes is 14.9 years


#----------------------------------------------------------------------------------------
#Figure 1: Kaplan-Meier curves
#----------------------------------------------------------------------------------------
pdf('Figure1.pdf', height = 11/2, width = 8.5/1.5)
  library(rms)
  km_rev <- npsurv(Surv(cus_t_surv, cus_fail) ~ cus_sr, data = cus_spells)
  survplot(km_rev, xlab = 'Years', xlim = c(0, 60), lty=c(2, 1), label.curves = FALSE, conf.int = .95)
  legend(0, 0.201, c('Revolutionary Regimes', 'Non-Revolutionary Regimes'),
         lty=c(1, 2),
         bty = "n",
         cex=0.7)
dev.off()

tiff(file = 'Figure1.tiff', height = 3*3200/2, width = 3*3200/2, units = 'px', res = 800)
  library(rms)
  km_rev <- npsurv(Surv(cus_t_surv, cus_fail) ~ cus_sr, data = cus_spells)
  survplot(km_rev, xlab = 'Years', xlim = c(0, 60), lty=c(2, 1), label.curves = FALSE, conf.int = .95)
  legend(0, 0.201, c('Revolutionary Regimes', 'Non-Revolutionary Regimes'),
         lty=c(1, 2),
         bty = "n",
         cex=0.7)
dev.off()
  
#Annual rates of breakdown 
sum(cus_spells$cus_fail[cus_spells$cus_sr==1])/
  sum(cus_spells$cus_t_surv[cus_spells$cus_sr==1])*100#1.28% for revolutionary regimes
sum(cus_spells$cus_fail[cus_spells$cus_sr==0])/
  sum(cus_spells$cus_t_surv[cus_spells$cus_sr==0])*100#5.80% for nonrevolutionary regimes

sum(cus_spells$cus_fail[cus_spells$cus_sr==1&cus_spells$partB==1])/
  sum(cus_spells$cus_t_surv[cus_spells$cus_sr==1&cus_spells$partB==1])*100#1.35% for revolutionary party regimes
sum(cus_spells$cus_fail[cus_spells$cus_sr==0&cus_spells$partB==1])/
  sum(cus_spells$cus_t_surv[cus_spells$cus_sr==0&cus_spells$partB==1])*100#2.95% for non-revolutionary party regimes

sum(cus_spells$cus_fail[cus_spells$cus_sr==1&cus_spells$sv_comm==1])/
  sum(cus_spells$cus_t_surv[cus_spells$cus_sr==1&cus_spells$sv_comm==1])*100#0.92% for revolutionary communist regimes
sum(cus_spells$cus_fail[cus_spells$cus_sr==0&cus_spells$sv_comm==1])/
  sum(cus_spells$cus_t_surv[cus_spells$cus_sr==0&cus_spells$sv_comm==1])*100#2.20% for non-revolutionary communist regimes

#Log rank tests 
diff.out <- survdiff(Surv(cus_t_surv, cus_fail) ~ cus_sr, data =  cus_spells)
pchisq(diff.out$chisq, length(diff.out$n)-1, lower.tail = FALSE)# p < 0.001 for rev VS non-rev

diff.out <- survdiff(Surv(cus_t_surv, cus_fail) ~ cus_sr, data = cus_spells, subset = partB == 1)
pchisq(diff.out$chisq, length(diff.out$n)-1, lower.tail = FALSE)# p = 0.0052 for rev party VS non-rev party

diff.out <- survdiff(Surv(cus_t_surv, cus_fail) ~ cus_sr, data = cus_spells, subset = sv_comm == 1)
pchisq(diff.out$chisq, length(diff.out$n)-1, lower.tail = FALSE)# p = 0.018 for rev communist VS non-rev communist 

#30-year survival rates 
km <- survfit(Surv(cus_t_surv, cus_fail) ~ 1, data = cus_spells, subset = cus_sr == 1)
  survest <- stepfun(km$time, c(1, km$surv))
  survest(30)#71% for revolutionary regimes 
km <- survfit(Surv(cus_t_surv, cus_fail) ~ 1, data = cus_spells, subset = cus_sr == 0)
  survest <- stepfun(km$time, c(1, km$surv))
  survest(30)#19% for revolutionary regimes 
  
#----------------------------------------------------------------------------------------
#Table 2: Summary Statistics
#----------------------------------------------------------------------------------------
varbles <- c('cus_sr', 'e_migdppclnL2', 'e_migdpgroL2', 'ln.tpopL2', 'ln.oil_gas_valuePOP_2014L2', 'partB', 'persB', 'mil', 'mon', 'sv_comm', 'InstalledCom')
sumStatsSpells <- stargazer(cus_spells[varbles], flip = FALSE, title = 'Summary Statistics, Fixed Covariates', notes.append = FALSE, align = TRUE, float = FALSE, digits = 2)
for (entry in dictSummstats){
  sumStatsSpells <- gsub(entry, names(dictSummstats)[dictSummstats == entry], sumStatsSpells)
}

  
#----------------------------------------------------------------------------------------
#Table 3: Cox regression
#----------------------------------------------------------------------------------------
form1 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2
model1.base <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = update(form1, . ~ . -cus_sr))
model1 <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = form1)
model1$summary; bet1 <- model1$summary[1,1]; se1 <- model1$summary[1,2];bet1
1-exp(bet1)#72 percent reduction in the hazards
1-exp(bet1-1.96*se1)#Upper CI: 85 percent reduction in the hazards
1-exp(bet1+1.96*se1)#Lower CI: 49 percent reduction in the hazards
model1$N
2*pnorm(-abs(bet1/se1))

form2 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 + ln.oil_gas_valuePOP_2014L2
model2.base <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = update(form2, . ~ . -cus_sr))
model2 <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = form2)
model2$summary; bet2 <- model2$summary[1,1]; se2 <- model2$summary[1,2];bet2
1-exp(bet2)
1-exp(bet2+1.96*se2)
1-exp(bet2-1.96*se2)
model2$N
2*pnorm(-abs(bet2/se2))

form3 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2+ ln.tpopL2 + ln.oil_gas_valuePOP_2014L2  + strata(e_regionpol)
model3.base <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = update(form3, . ~ . -cus_sr))
model3 <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = form3)
model3$summary; bet3 <- model3$summary[1,1]; se3 <- model3$summary[1,2];bet3
1-exp(bet3)
1-exp(bet3+1.96*se3)
1-exp(bet3-1.96*se3)
model3$N
2*pnorm(-abs(bet3/se3))

form4 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 + ln.oil_gas_valuePOP_2014L2 + strata(cus_ARevnum)
model4.base <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = update(form4, . ~ . -cus_sr))
model4 <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = form4)
model4$summary; bet4 <- model4$summary[1,1]; se4 <- model4$summary[1,2]; bet4
1-exp(bet4)# decrease in the hazards of 80%
1-exp(bet4+1.96*se4)
1-exp(bet4-1.96*se4)
model4$N
2*pnorm(-abs(bet4/se4))

form5 <- Surv(cus_t_surv, cus_fail) ~ cus_sr+ e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 + ln.oil_gas_valuePOP_2014L2 + partB +persB + mil + mon
model5.base <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = update(form5, . ~ . -cus_sr))
model5 <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = form5)
model5$summary; bet5 <- model5$summary[1,1]; se5 <- model5$summary[1,2];bet5
1-exp(bet5)#origins are associated with a 58 percent decrease in the hazards in Model 5
1-exp(bet5+1.96*se5)
1-exp(bet5-1.96*se5)
model5$N
2*pnorm(-abs(bet5/se5))

#Party and Monarchy are associated with a decrease of 64% and 81%
1-exp(-1.0179)#part
1-exp(-1.661)#mon

form6 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 + ln.oil_gas_valuePOP_2014L2 + partB +persB + mil + mon  + strata(cus_ARevnum)
model6.base <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = update(form6, . ~ . -cus_sr))
model6 <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = form6)
model6$summary; bet6 <- model6$summary[1,1]; se6 <- model6$summary[1,2]; bet6
1-exp(bet6)#origins are associated with a 67 percent decrease in the hazards in Model 6
1-exp(bet6+1.96*se6)
1-exp(bet6-1.96*se6)
model6$N
2*pnorm(-abs(bet6/se6))

form7 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2  + ln.oil_gas_valuePOP_2014L2+ partB +persB + mil + mon + sv_comm + strata(cus_ARevnum)
model7.base <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = update(form7, . ~ . -cus_sr))
model7 <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = form7)
model7$summary; bet7 <- model7$summary[1,1]; se7 <- model7$summary[1,2]
1-exp(bet7)#Revolution is associated with a 62 percent decrease in the hazards in Model 7 (p < 0.05)
1-exp(bet7+1.96*se7)
1-exp(bet7-1.96*se7)
model7$N
2*pnorm(-abs(bet7/se7))

form8 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 +ln.oil_gas_valuePOP_2014L2  +  partB +persB + mil + mon + sv_comm+ InstalledCom + strata(cus_ARevnum)
model8.base <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = update(form8, . ~ . -cus_sr))
model8 <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = form8)
model8$summary; bet8 <- model8$summary[1,1]; se8 <- model8$summary[1,2]
1-exp(bet8)#Revolution is associated with a 70 percent decrease in the hazards when we also condition on the foreign imposition of communism.
1-exp(bet8+1.96*se8)
1-exp(bet8-1.96*se8)
model8$N
2*pnorm(-abs(bet8/se8))

ref.models <- list(#models without controlling for revolution
  model1.base,
  model2.base,
  model3.base,
  model4.base,
  model5.base,
  model6.base,
  model7.base,
  model8.base
)

modelListObject <-
  CoxmodelList(
    model1,
    model2,
    model3,
    model4,
    model5,
    model6,
    model7,
    model8
  )

ModelLabels <- c('Non-stratified',
                 'Stratified',
                 'Controlling for regime type (post-treatment)',
                 'Controlling for communism (post-treatment)')
#library(plyr)
Table3 <- CoxTable(modelListObject = modelListObject,
                      ref.Models = ref.models,
                      dict = dictSummstats,
                      column.labels = ModelLabels,
                      column.separate = c(2, 2, 2, 2),
                      dep.var.labels = 'Authoritarian Regime Breakdown',
                      align = TRUE,
                      float = FALSE,
                      font.size = 'footnotesize')


Table3_htlm <- CoxTable(modelListObject = modelListObject,
                   ref.Models = ref.models,
                   dict = dictSummstats,
                   column.labels = ModelLabels,
                   column.separate = c(2, 2, 2, 2),
                   dep.var.labels = 'Authoritarian Regime Breakdown',
                   align = TRUE,
                   float = FALSE,
                   font.size = 'footnotesize',
                   type.in = 'html')

#LR tests
#coxLRtest1 <- coxLRtest(fit1 = modelListObject$models[[1]], fit0 = ref.models[[1]])# p = 0.000048
#coxLRtest2 <- coxLRtest(fit1 = modelListObject$models[[2]], fit0 = ref.models[[2]])# p = 0.00014
#coxLRtest3 <- coxLRtest(fit1 = modelListObject$models[[3]], fit0 = ref.models[[3]])# p = 0.0000048
#coxLRtest4 <- coxLRtest(fit1 = modelListObject$models[[4]], fit0 = ref.models[[4]])# p = 0.000008
#coxLRtest5 <- coxLRtest(fit1 = modelListObject$models[[5]], fit0 = ref.models[[5]])# p = 0.013
#coxLRtest6 <- coxLRtest(fit1 = modelListObject$models[[6]], fit0 = ref.models[[6]])# p = 0.007
#coxLRtest7 <- coxLRtest(fit1 = modelListObject$models[[7]], fit0 = ref.models[[7]])# p = 0.027
#coxLRtest8 <- coxLRtest(fit1 = modelListObject$models[[8]], fit0 = ref.models[[8]])# p = 0.012

#footnote 121
#Analysis of deviance
cox.1 <- coxph( Surv(cus_t_surv, cus_fail) ~ cus_sr + partB + persB + mil + mon, data = cus_spells, robust=T)
cox.0 <- coxph( Surv(cus_t_surv, cus_fail) ~ partB +persB + mil + mon, data = cus_spells, robust=T)
anova(cox.1, cox.0, test = 'Chisq')

#----------------------------------------------------------------------------------------
#Figure 2
#----------------------------------------------------------------------------------------
form_gwf <- formula(Surv(gwf_t_surv, gwf_fail) ~ gwf_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 + ln.oil_gas_valuePOP_2014L2  )
model_gwf <- ameliaCoxEstimate(ameliaObject = d_mi_spells_gwf, formul = form_gwf, fail_var = 'gwf_fail', id_var = 'gwf_caseid')
model_gwf$summary
model_gwf$N
summary(coxph(form_gwf, d_mi_spells_gwf$imputations[[1]]), robust=T)
sum(d_mi_spells_gwf$imputations[[1]]$gwf_sr)

#For Svolik, only one lag is necessary
form_svol <- formula(Surv(svol_t_surv, svol_fail) ~ svol_sr + e_migdppclnL1 + e_migdpgroL1 + ln.tpopL1 + ln.oil_gas_valuePOP_2014L1)
model_svol <- ameliaCoxEstimate(ameliaObject = d_mi_spells_svol, formul = form_svol, fail_var = 'svol_fail', id_var = 'svol_caseid')
model_svol$summary
model_svol$N
summary(coxph(form_svol, d_mi_spells_svol$imputations[[1]]), robust=T)
sum(d_mi_spells_svol$imputations[[1]]$svol_sr)

form_mcm <- formula(Surv(mcm_t_surv, mcm_fail) ~ mcm_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 + ln.oil_gas_valuePOP_2014L2 )
model_mcm <- ameliaCoxEstimate(ameliaObject = d_mi_spells_mcm, formul = form_mcm, fail_var = 'mcm_fail', id_var = 'mcm_caseid')
model_mcm$summary
model_mcm$N
summary(coxph(form_mcm, d_mi_spells_mcm$imputations[[1]]), robust=T)
sum(d_mi_spells_mcm$imputations[[1]]$mcm_sr)

form_N20 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 + ln.oil_gas_valuePOP_2014L2
model_N20 <- ameliaCoxEstimate(ameliaObject = d_mi_spells_N20, formul = form_N20, id_var = 'cus_caseid')
model_N20$summary
model_N20$N
summary(coxph(form_N20, d_mi_spells_N20$imputations[[1]]), robust=T)
sum(d_mi_spells_N20$imputations[[1]]$cus_sr)

coefs <- c(
  model_gwf$beta[1],
  model_svol$beta[1],
  model_mcm$beta[1],
  model_N20$beta[1]
)

ses <- c(
  model_gwf$se[1],
  model_svol$se[1],
  model_mcm$se[1],
  model_N20$se[1]
)

the_names <- c('GWF', 'Svolik', 'MCM', 'Expanded' )
lower95 <- qnorm(0.025, mean = coefs, sd = ses)
upper95 <- qnorm(0.975, mean = coefs, sd = ses)
the_dat <- data.frame(coefs = coefs,
                      lower95 = lower95,
                      upper95 = upper95
)

pdf('Figure2.pdf', width = '7in', height = '3in')
  plot(x = seq(0, length(the_names) - 1), y = coefs,
       ylab = 'Revolution coefficient',
       xlab= '',
       xaxt = 'n',
       yaxt = 'n',
       xlim = c(-1, length(the_names)),
       ylim = c(-2.5, 1.3),
       las= 1)
  abline(0, 0, lty = 3)
  axis(side =1, line = 1, at = 0:(length(the_names) - 1), labels = the_names, cex.lab = 0.5)
  axis(las= 1, side =2, at = -3:3, labels = -3:3, cex.lab = 0.5)
  error_bars <- arrows(x0 = 0:(length(the_names) - 1), y0 = lower95, x1 = 0:(length(the_names) - 1), y1 = upper95, code = 3, angle = 90, length=0.05)
dev.off()

tiff(file = 'Figure2.tiff', height = 3*3200/2, width = 3*3200/2, units = 'px', res = 800)
  plot(x = seq(0, length(the_names) - 1), y = coefs,
       ylab = 'Revolution coefficient',
       xlab= '',
       xaxt = 'n',
       yaxt = 'n',
       xlim = c(-1, length(the_names)),
       ylim = c(-2.5, 1.3),
       las= 1)
  abline(0, 0, lty = 3)
  axis(side =1, line = 1, at = 0:(length(the_names) - 1), labels = the_names, cex.lab = 0.5)
  axis(las= 1, side =2, at = -3:3, labels = -3:3, cex.lab = 0.5)
  error_bars <- arrows(x0 = 0:(length(the_names) - 1), y0 = lower95, x1 = 0:(length(the_names) - 1), y1 = upper95, code = 3, angle = 90, length=0.05)
dev.off()

#----------------------------------------------------------------------------------------
#CBPS Reweighing
#----------------------------------------------------------------------------------------
library(CBPS)
library(ggplot2)
library(cobalt)
mi.vars <- c(
  'e_migdppclnL2',
  'e_migdpgroL2',
  'warL2',
  'prev_demo',
  'prev_partB',
  'prev_persB',
  'prev_mil',
  'prev_mon'
)


match.frml <- as.formula(paste('cus_sr ~ ', paste(mi.vars, collapse = '+'), sep = ''))
dd.mi <- NULL
for (i in 1:d_mi_spells_cus$m) {
  dd <- d_mi_spells_cus$imputations[[i]]
  dd$cus_fail <- cus_spells$cus_fail
  dd$cus_t_surv <- cus_spells$cus_t_surv
  dd$Country <- cus_spells$Country
  dd$prev_demo <- cus_spells$prev_demo
  dd$prev_mil <- cus_spells$prev_mil
  dd$prev_partB <- cus_spells$prev_partB
  dd$prev_mon <- cus_spells$prev_mon
  dd$prev_persB <- cus_spells$prev_persB
  m.out.att <- CBPS(match.frml, ATT=1, data = dd)
  dd$w <- m.out.att$weights
  dd$.imp <- i
  dd.mi <- rbind(dd.mi, dd)
}

library(survey)
cfs1 <- NULL
ses1 <- NULL
cfs2 <- NULL
ses2 <- NULL
for (i in 1:d_mi_spells_cus$m){
  the_d <- subset(dd.mi, .imp==i)
  the_d$id <- as.character(1:nrow(the_d))
  des <- svydesign(ids = ~1, weights = the_d$w,
                   data = the_d)
  svycoxph.out <- svycoxph(Surv(cus_t_surv, cus_fail)~cus_sr +cluster(id),design=des)
  cox.out <- coxph(Surv(cus_t_surv, cus_fail)~cus_sr+cluster(id), data = the_d, weights = the_d$w)
  summary(cox.out)
  cfs1 <- rbind(cfs1, summary(svycoxph.out)$coefficients[1, 1])
  ses1 <- rbind(ses1, sqrt(diag(svycoxph.out$var))[1])
  cfs2 <- rbind(cfs2, summary(cox.out)$coefficients[1, 1])
  ses2 <- rbind(ses2, sqrt(diag(cox.out$var))[1])
}

beta1 <- mi.meld(cfs1, ses1)$q.mi; se1 <- mi.meld(cfs1, ses1)$se.mi;beta1
z1 <- beta1/se1
2*(pnorm(-abs(z1)))
beta2 <- mi.meld(cfs2, ses2)$q.mi; se1 <- mi.meld(cfs2, ses2)$se.mi;beta2
z2 <- beta2/se2
2*pnorm(-abs(z2))
1-exp(beta1)
1-exp(beta1+1.96*se1)
1-exp(beta1-1.96*se1)

#Covariate balance plot
newNames <- c(
  'e_migdppclnL2' = 'GDP per capita (log)',
  'e_migdpgroL2' = 'GDP growth',
  'ln.tpopL2' = 'Population size (log)',
  'warL2' = 'War',
  'prev_demo' = 'Previous democracy',
  'prev_partB' = 'Previous party',
  'prev_persB' = 'Previous personalist',
  'prev_mil' = 'Previous military',
  'prev_mon' = 'Previous monarchy'
)

balances <- bal.tab(match.frml,
                    data = dd.mi, weights = "w", imp = ".imp", method = "weighting")

library(ggplot2); library(cobalt)
figA4 <- cobalt::love.plot(balances,
          agg.fun = 'range',
          threshold = .25,
          var.order = 'unadjusted',
          stat = 'mean.diffs',
          line = F,
          abs = F,
          limits = c(-0.65, 0.65),
          #limits = c(0, 3),
          var.names = newNames
          #shapes = c(22, 24),
) +
  geom_vline(xintercept=.1, linetype = 2, color = "black") +
  geom_vline(xintercept=-.1, linetype = 2, color = "black") #+
ggsave(filename = 'FigureA4.pdf', plot = figA4)

#-------------------------------
#Coarsened exact matching
#-------------------------------
library(cem)
d.cem <- cus_spells

m.out <- cem('cus_sr',
             data=cus_spells[, c('cus_sr', 'cus_t_surv', 'cus_fail', mi.vars)],
             drop = c('cus_sr', 'cus_t_surv', 'cus_fail', 'partB',
            baseline = 0,
              eval.imbalance = T)
)
cus_spells_matched <- cus_spells[m.out$matched, ]; nrow(cus_spells_matched)#70 regimes matched
sum(cus_spells_matched$cus_sr)#12 revolutionary regimes
sum(!cus_spells_matched$cus_sr)#58 non-revolutionary regimes
cus_spells_matched$strata <- m.out$strata[m.out$matched]#strata
cus_spells_matched$w <- m.out$w[m.out$matched]#weights
cus_spells_matched$id <- as.factor(1:nrow(cus_spells_matched))#het-consistent standard errors

fsatt.out <- coxph(Surv(cus_t_surv, cus_fail)~cus_sr+cluster(id),data=cus_spells_matched, weights = w)
beta <- summary(fsatt.out)$coefficients[1,1];beta
se <- summary(fsatt.out)$coefficients[1,4];se
2*pnorm(-abs(beta/se))

#Footnote 138 -> list revolutionary cases from matching
as.character(cus_spells_matched$cus_caseid[cus_spells_matched$cus_sr==1])

####################################################################
#Table 4: Mediators
####################################################################
cus_spells$sdate <- apply(cus_spells[, c('cus_syear', 'cus_smonth', 'cus_sday')], 1, function(x) paste(x, collapse = '-'))
cus_spells$edate <- apply(cus_spells[, c('cus_eyear', 'cus_emonth', 'cus_eday')], 1, function(x) paste(x, collapse = '-'))
cus_spells$edate[cus_spells$edate == 'NA-NA-NA'] <- ''

#Civil society strength
library(ri); set.seed(93723)
d.civsoc <- na.omit(cus_spells[, c('cus_sr', 'ccode', 'v2xcs_ccsi_avg')])
civsoc.p.val <- r_inference(d.civsoc, 'v2xcs_ccsi_avg')$two.tailed.p.value
civsoc.p.val
row.civsoc <- c(
  'Civil society strength',
  deci(mean(d.civsoc$v2xcs_ccsi_avg[d.civsoc$cus_sr == 0], na.rm=T), 3),
  deci(mean(d.civsoc$v2xcs_ccsi_avg[d.civsoc$cus_sr == 1], na.rm=T), 3),
  deci(mean(d.civsoc$v2xcs_ccsi_avg[d.civsoc$cus_sr == 1], na.rm=T) - mean(d.civsoc$v2xcs_ccsi_avg[d.civsoc$cus_sr == 0], na.rm = T), 3),
  paste(deci(civsoc.p.val, 3), sep = '')
)

#Military size
cus_spells$ln.milperpc_avg <- log(cus_spells$milperpc_avg)
d.milper <- na.omit(cus_spells[is.finite(cus_spells$ln.milperpc_avg), c('cus_sr', 'ccode', 'ln.milperpc_avg')])
d.milper.p.val <- r_inference(some_data = d.milper, dv = 'ln.milperpc_avg')$two.tailed.p.value
row.milper <- c(
  'Military size per capita (log)',
  deci(mean(d.milper$ln.milperpc_avg[d.milper$cus_sr == 0], na.rm=T), 3),
  deci(mean(d.milper$ln.milperpc_avg[d.milper$cus_sr == 1], na.rm=T), 3),
  deci(mean(d.milper$ln.milperpc_avg[d.milper$cus_sr == 1], na.rm=T) - mean(d.milper$ln.milperpc_avg[d.milper$cus_sr == 0], na.rm = T), 3),
  paste(deci(d.milper.p.val, 3), sep = '')
)

#Party milit 
d.partymilit <- na.omit(cus_spells[, c('cus_sr', 'ccode', 'partymilit_avg')])
partymilit.p.val <- r_inference(d.partymilit, 'partymilit_avg')$two.tailed.p.value
row.partymilit <- c(
  'Party control over the military',
  deci(mean(d.partymilit$partymilit_avg[d.partymilit$cus_sr == 0], na.rm=T), 3),
  deci(mean(d.partymilit$partymilit_avg[d.partymilit$cus_sr == 1], na.rm=T), 3),
  deci(mean(d.partymilit$partymilit_avg[d.partymilit$cus_sr == 1], na.rm=T) - mean(d.partymilit$partymilit_avg[d.partymilit$cus_sr == 0], na.rm = T), 3),
  '$<0.001$'
)

#Coups
d.coups <- na.omit(cus_spells[, c('cus_sr', 'ccode', 'pth_att_avg')])
pth.att.p.val <- r_inference(d.coups, 'pth_att_avg')$two.tailed.p.value
row.cp.pth <- c(
  'Annual rate of coup attempts',
  deci(mean(d.coups$pth_att_avg[d.coups$cus_sr == 0], na.rm=T), 3),
  deci(mean(d.coups$pth_att_avg[d.coups$cus_sr == 1], na.rm=T), 3),
  deci(mean(d.coups$pth_att_avg[d.coups$cus_sr == 1], na.rm=T) - mean(d.coups$pth_att_avg[d.coups$cus_sr == 0], na.rm = T), 3),
  paste(deci(pth.att.p.val, 3), sep = '')
)

#Party rule
row.part <- c(
  'Presence of ruling party',
  deci(mean(cus_spells$partB[cus_spells$cus_sr == 0]), 3),
  deci(mean(cus_spells$partB[cus_spells$cus_sr == 1]), 3),
  deci(mean(cus_spells$partB[cus_spells$cus_sr == 1]) - mean(cus_spells$partB[cus_spells$cus_sr == 0]), 3),
  '$< 0.001$'
)
tab4 <- matrix("", nrow = 0, ncol = 4)
row0 <- c('Mediator', 'Non Rev.', 'Rev.', 'Difference', 'p-value')
row0a <- c('', '', '', '', 'on difference')
row0c <- c('', '', '', '', '')

tab4 <- rbind(
  row0,
  row.civsoc,
  row.milper,
  row.partymilit,
  row.cp.pth,
  row.part
)

library(xtable)
Table4 <- print(
  xtable(tab4,
         caption = "Mediator Variables, Differences in Means Tests",
         align = 'l|l|cc|cc|',
         label = 'balTab'),
  include.rownames = FALSE,
  caption.placement = 'top',
  table.placement = "H",
  include.colnames = FALSE,
  size = "\\scriptsize",
  hline.after = c(0, 1, 1, 2, 3, 5, 6),
  sanitize.text.function = function(x){x})
Table4 <- gsub('-', '\u2013', Table4)
Table4 <- gsub("\\@\\@",  "\\ \\ \\ \\", Table4)

#Footnote 147
#Median military size
median(cus_spells[cus_spells$cus_sr ==1, 'milperpc_avg'], na.rm=T)*1000
median(cus_spells[cus_spells$cus_sr ==0, 'milperpc_avg'], na.rm=T)*1000

####################################################################
#Table 5: Mechanism ATTs
####################################################################
match.frml <- as.formula(paste('cus_sr ~ ', paste(mi.vars, collapse = '+'), sep = ''))
d.out.civsoc <- NULL
d.out.milper <- NULL
d.out.partymilit <- NULL
d.out.coups <- NULL
d.out.party <- NULL
for (i in 1:d_mi_spells_cus$m) {
  d.mechs <- d_mi_spells_cus$imputations[[i]]
  d.mechs$prev_demo <- cus_spells$prev_demo
  d.mechs$prev_partB <- cus_spells$prev_partB
  d.mechs$prev_persB <- cus_spells$prev_persB
  d.mechs$prev_mil <- cus_spells$prev_mil
  d.mechs$prev_mon <- cus_spells$prev_mon
  
  d.mechs$v2xcs_ccsi_avg <- cus_spells$v2xcs_ccsi_avg
  d.mechs$ln.milperpc_avg <- cus_spells$ln.milperpc_avg
  d.mechs$partymilit_avg <- cus_spells$partymilit_avg
  d.mechs$pth_att_avg <- cus_spells$pth_att_avg

  d.mechs.civsoc <- na.omit(d.mechs[, c('v2xcs_ccsi_avg', all.vars(match.frml))])
  d.mechs.milper <- na.omit(d.mechs[, c('ln.milperpc_avg', all.vars(match.frml))]);d.mechs.milper <- d.mechs.milper[is.finite(d.mechs.milper$ln.milperpc_avg), ]
  d.mechs.partymilit <- na.omit(d.mechs[, c('partymilit_avg', all.vars(match.frml))])
  d.mechs.coups <- na.omit(d.mechs[, c('pth_att_avg', all.vars(match.frml))])
  d.mechs.party <- na.omit(d.mechs[, c('partB', all.vars(match.frml))])
  
  cbps.out.civsoc <- CBPS(match.frml, ATT = 1, data = d.mechs.civsoc)
  cbps.out.milper <- CBPS(match.frml, ATT = 1, data = d.mechs.milper)
  cbps.out.partymilit <- CBPS(match.frml, ATT = 1, data = d.mechs.partymilit)
  cbps.out.coups <- CBPS(match.frml, ATT = 1, data = d.mechs.coups)
  cbps.out.party <- CBPS(match.frml, ATT = 1, data = d.mechs.party)
  
  d.mechs.civsoc$w <- cbps.out.civsoc$weight
  d.mechs.milper$w <- cbps.out.milper$weight
  d.mechs.partymilit$w <- cbps.out.partymilit$weight
  d.mechs.coups$w <- cbps.out.coups$weight
  d.mechs.party$w <- cbps.out.party$weight
  
  d.mechs.civsoc$.imp <- i
  d.mechs.milper$.imp <- i
  d.mechs.partymilit$.imp <- i
  d.mechs.coups$.imp <- i
  d.mechs.party$.imp <- i
  
  d.out.civsoc <- rbind(d.out.civsoc, d.mechs.civsoc)
  d.out.milper <- rbind(d.out.milper, d.mechs.milper)
  d.out.partymilit <- rbind(d.out.partymilit, d.mechs.partymilit)
  d.out.coups <- rbind(d.out.coups, d.mechs.coups)
  d.out.party <- rbind(d.out.party, d.mechs.party)
}

library(lmtest)
library(sandwich)
cfs.civsoc <- ses.civsoc <- NULL
cfs.milper <- ses.milper<- NULL
cfs.partymilit <- ses.partymilit<- NULL
cfs.coup <- ses.coup<- NULL
cfs.party <- ses.party<- NULL
for (i in 1:d_mi_spells_cus$m){
  d.civsoc.att <- subset(d.out.civsoc, .imp==i)
  d.milper.att <- subset(d.out.milper, .imp==i)
  d.partymilit.att <- subset(d.out.partymilit, .imp==i)
  d.coup.att <- subset(d.out.coups, .imp==i)
  d.party.att <- subset(d.out.party, .imp==i)
  
  lm.civsoc <- lm(v2xcs_ccsi_avg~cus_sr , data = d.civsoc.att, weights = w)
  lm.milper <- lm(ln.milperpc_avg~cus_sr , data = d.milper.att, weights = w)
  lm.partymilit <- lm(partymilit_avg~cus_sr , data = d.partymilit.att, weights = w)
  lm.coup <- lm(pth_att_avg~cus_sr , data = d.coup.att, weights = w)
  lm.party <- lm(partB~cus_sr , data = d.party.att, weights = w)
  
  cfs.civsoc <- rbind(cfs.civsoc, coeftest(lm.civsoc, df = df.residual(lm.civsoc), vcov = vcovHC, type = "HC3")[2,1])
  ses.civsoc <- rbind(ses.civsoc, coeftest(lm.civsoc, df = df.residual(lm.civsoc), vcov = vcovHC, type = "HC3")[2,2])
  
  cfs.milper <- rbind(cfs.milper, coeftest(lm.milper, df = df.residual(lm.milper), vcov = vcovHC, type = "HC3")[2,1])
  ses.milper <- rbind(ses.milper, coeftest(lm.milper, df = df.residual(lm.milper), vcov = vcovHC, type = "HC3")[2,2])
  
  cfs.partymilit <- rbind(cfs.partymilit, coeftest(lm.partymilit, df = df.residual(lm.partymilit), vcov = vcovHC, type = "HC3")[2,1])
  ses.partymilit <- rbind(ses.partymilit, coeftest(lm.partymilit, df = df.residual(lm.partymilit), vcov = vcovHC, type = "HC3")[2,2])
  
  cfs.coup <- rbind(cfs.coup, coeftest(lm.coup, df = df.residual(lm.coup), vcov = vcovHC, type = "HC3")[2,1])
  ses.coup <- rbind(ses.coup, coeftest(lm.coup, df = df.residual(lm.coup), vcov = vcovHC, type = "HC3")[2,2])
  
  cfs.party <- rbind(cfs.party, coeftest(lm.party, df = df.residual(lm.party), vcov = vcovHC, type = "HC3")[2,1])
  ses.party <- rbind(ses.party, coeftest(lm.party, df = df.residual(lm.party), vcov = vcovHC, type = "HC3")[2,2])
}

civsoc.out <- mi.meld(cfs.civsoc, ses.civsoc)
milper.out <- mi.meld(cfs.milper, ses.milper)
partymilit.out <- mi.meld(cfs.partymilit, ses.partymilit)
coup.out <- mi.meld(cfs.coup, ses.coup)
party.out <- mi.meld(cfs.party, ses.party)

#
#This effect is substantively large and corresponds to about 0.4 standard deviation.
civsoc.out$q.mi/sd(dat$v2xcs_ccsi, na.rm=T)

mediators <- c(
  'Mediator',
  'Civil society strength',
  '',
  'Military size per capita (log)',
  '',
  'Party control over the military',
  '',
  'Annual rate of coup attempts',
  '',
  'Party rule',
  ''
)

atts <- c(
  'ATT',
  deci(civsoc.out$q.mi, 3),
  paste('(', deci(civsoc.out$se.mi, 3), ')', sep=''),
  deci(milper.out$q.mi, 3),
  paste('(', deci(milper.out$se.mi, 3), ')', sep=''),
  deci(partymilit.out$q.mi, 3),
  paste('(', deci(partymilit.out$se.mi, 3), ')', sep=''),
  deci(coup.out$q.mi, 3),
  paste('(', deci(coup.out$se.mi, 3), ')', sep=''),
  deci(party.out$q.mi, 3),
  paste('(', deci(party.out$se.mi, 3), ')', sep='')
)

#format(p_vals)
options(scipen=999)
p_vals <- c(
  2*pnorm(-abs(civsoc.out$q.mi/civsoc.out$se.mi)),
  2*pnorm(-abs(milper.out$q.mi/milper.out$se.mi)),
  2*pnorm(-abs(partymilit.out$q.mi/partymilit.out$se.mi)),
  2*pnorm(-abs(coup.out$q.mi/coup.out$se.mi)),
  2*pnorm(-abs(party.out$q.mi/party.out$se.mi))
)

names(p_vals) <- c('civsoc', 'milper', 'partymilit', 'coup', 'party')

exp(.933)

library(xtable)
Table5 <- print(
  xtable(xtable(cbind(mediators, atts)),
         caption = "Effect of Revolution on the Mediator Variables (ATT)",
         align = 'l|l|c|',
         label = 'balTab'),
  include.rownames = FALSE,
  caption.placement = 'top',
  table.placement = "H",
  include.colnames = FALSE,
  size = "\\scriptsize",
  hline.after = c(0, 1, 1, 3, 5, 9, 11),
  sanitize.text.function = function(x){x})

Table5 <- gsub('-', '\u2013', Table5)
Table5 <- gsub("\\@\\@",  "\\ \\ \\ \\", Table5)

####################################################################
#Figure 3: Synthetic control analysis
####################################################################

load('synth.Rdata')
pdf('Figure3.pdf', height = 8*2, width = 6*2)
  par(mfrow=c(5, 4), oma = c(0.001, 0.001, 0.001, 0.001))#
  title_af = paste('AFGHANISTAN', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Afghanistan, data.civsoc.Afghanistan, 4), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Afghanistan, data.civsoc.Afghanistan), 3), sep = '')
  data.civsoc.Afghanistan$Y1plot[data.civsoc.Afghanistan$tag$time.plot >4] <- NA
  data.civsoc.Afghanistan$Y0plot[data.civsoc.Afghanistan$tag$time.plot >4] <- NA
  path.plot.cus(dataprep.res = data.civsoc.Afghanistan, synth.res = synth.civsoc.Afghanistan, tr.intake = 1996, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)

    
  title_af = paste('ALBANIA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Albania, data.civsoc.Albania), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Albania, data.civsoc.Albania), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Albania, synth.res = synth.civsoc.Albania, tr.intake = 1944, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('ALGERIA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Algeria, data.civsoc.Algeria), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Algeria, data.civsoc.Algeria), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Algeria, synth.res = synth.civsoc.Algeria, tr.intake = 1962, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)

  title_af = paste('ANGOLA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Angola, data.civsoc.Angola), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Angola, data.civsoc.Angola), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Angola, synth.res = synth.civsoc.Angola, tr.intake = 1975, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('BOLIVIA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Bolivia, data.civsoc.Bolivia), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Bolivia, data.civsoc.Bolivia), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Bolivia, synth.res = synth.civsoc.Bolivia, tr.intake = 1952, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('CAMBODIA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Cambodia, data.civsoc.Cambodia,3), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Cambodia, data.civsoc.Cambodia), 3), sep = '')
  data.civsoc.Cambodia$Y1plot[data.civsoc.Cambodia$tag$time.plot >3] <- NA
  data.civsoc.Cambodia$Y0plot[data.civsoc.Cambodia$tag$time.plot >3] <- NA
  path.plot.cus(dataprep.res = data.civsoc.Cambodia, synth.res = synth.civsoc.Cambodia, tr.intake = 1975, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('CHINA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.China, data.civsoc.China), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.China, data.civsoc.China), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.China, synth.res = synth.civsoc.China, tr.intake = 1949, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('CUBA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Cuba, data.civsoc.Cuba), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Cuba, data.civsoc.Cuba), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Cuba, synth.res = synth.civsoc.Cuba, tr.intake = 1959, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('ERITREA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Eritrea, data.civsoc.Eritrea), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Eritrea, data.civsoc.Eritrea), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Eritrea, synth.res = synth.civsoc.Eritrea, tr.intake = 1993, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('GUINEA-BISSAU', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.GuineaBissau, data.civsoc.GuineaBissau), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.GuineaBissau, data.civsoc.GuineaBissau), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.GuineaBissau, synth.res = synth.civsoc.GuineaBissau, tr.intake = 1974, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('IRAN', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Iran, data.civsoc.Iran), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Iran, data.civsoc.Iran), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Iran, synth.res = synth.civsoc.Iran, tr.intake = 1979, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('MEXICO', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Mexico, data.civsoc.Mexico), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Mexico, data.civsoc.Mexico), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Mexico, synth.res = synth.civsoc.Mexico, tr.intake = 1915, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  
  title_af = paste('MOZAMBIQUE', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Mozambique, data.civsoc.Mozambique), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Mozambique, data.civsoc.Mozambique), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Mozambique, synth.res = synth.civsoc.Mozambique, tr.intake = 1975, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('NICARAGUA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Nicaragua, data.civsoc.Nicaragua), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Nicaragua, data.civsoc.Nicaragua), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Nicaragua, synth.res = synth.civsoc.Nicaragua, tr.intake = 1979, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('RWANDA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Rwanda, data.civsoc.Rwanda), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Rwanda, data.civsoc.Rwanda), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Rwanda, synth.res = synth.civsoc.Rwanda, tr.intake = 1994, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('U.S.S.R.', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Russia, data.civsoc.Russia), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Russia, data.civsoc.Russia), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Russia, synth.res = synth.civsoc.Russia, tr.intake = 1917, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('VIETNAM', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Vietnam, data.civsoc.Vietnam), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Vietnam, data.civsoc.Vietnam), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Vietnam, synth.res = synth.civsoc.Vietnam, tr.intake = 1954, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('YUGOSLAVIA 1945', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Yugoslavia, data.civsoc.Yugoslavia), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Yugoslavia, data.civsoc.Yugoslavia), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Yugoslavia, synth.res = synth.civsoc.Yugoslavia, tr.intake = 1945, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
dev.off()

#save to tiff

tiff(file = 'Figure3.tiff', height = 8*3200/2, width = 6*3200/2, units = 'px', res = 800)
  par(mfrow=c(5, 4), oma = c(0.001, 0.001, 0.001, 0.001))#
  #par(mar=c(1,1,1,1))
  title_af = paste('AFGHANISTAN', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Afghanistan, data.civsoc.Afghanistan, 4), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Afghanistan, data.civsoc.Afghanistan), 3), sep = '')
  data.civsoc.Afghanistan$Y1plot[data.civsoc.Afghanistan$tag$time.plot >4] <- NA
  data.civsoc.Afghanistan$Y0plot[data.civsoc.Afghanistan$tag$time.plot >4] <- NA
  path.plot.cus(dataprep.res = data.civsoc.Afghanistan, synth.res = synth.civsoc.Afghanistan, tr.intake = 1996, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  
  title_af = paste('ALBANIA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Albania, data.civsoc.Albania), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Albania, data.civsoc.Albania), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Albania, synth.res = synth.civsoc.Albania, tr.intake = 1944, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('ALGERIA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Algeria, data.civsoc.Algeria), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Algeria, data.civsoc.Algeria), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Algeria, synth.res = synth.civsoc.Algeria, tr.intake = 1962, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  data.civsoc.Angola$tag$time.plot
  
  title_af = paste('ANGOLA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Angola, data.civsoc.Angola), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Angola, data.civsoc.Angola), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Angola, synth.res = synth.civsoc.Angola, tr.intake = 1975, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('BOLIVIA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Bolivia, data.civsoc.Bolivia), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Bolivia, data.civsoc.Bolivia), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Bolivia, synth.res = synth.civsoc.Bolivia, tr.intake = 1952, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('CAMBODIA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Cambodia, data.civsoc.Cambodia,3), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Cambodia, data.civsoc.Cambodia), 3), sep = '')
  data.civsoc.Cambodia$Y1plot[data.civsoc.Cambodia$tag$time.plot >3] <- NA
  data.civsoc.Cambodia$Y0plot[data.civsoc.Cambodia$tag$time.plot >3] <- NA
  path.plot.cus(dataprep.res = data.civsoc.Cambodia, synth.res = synth.civsoc.Cambodia, tr.intake = 1975, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('CHINA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.China, data.civsoc.China), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.China, data.civsoc.China), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.China, synth.res = synth.civsoc.China, tr.intake = 1949, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('CUBA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Cuba, data.civsoc.Cuba), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Cuba, data.civsoc.Cuba), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Cuba, synth.res = synth.civsoc.Cuba, tr.intake = 1959, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('ERITREA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Eritrea, data.civsoc.Eritrea), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Eritrea, data.civsoc.Eritrea), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Eritrea, synth.res = synth.civsoc.Eritrea, tr.intake = 1993, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('GUINEA-BISSAU', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.GuineaBissau, data.civsoc.GuineaBissau), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.GuineaBissau, data.civsoc.GuineaBissau), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.GuineaBissau, synth.res = synth.civsoc.GuineaBissau, tr.intake = 1974, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('IRAN', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Iran, data.civsoc.Iran), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Iran, data.civsoc.Iran), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Iran, synth.res = synth.civsoc.Iran, tr.intake = 1979, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('MEXICO', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Mexico, data.civsoc.Mexico), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Mexico, data.civsoc.Mexico), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Mexico, synth.res = synth.civsoc.Mexico, tr.intake = 1915, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  
  title_af = paste('MOZAMBIQUE', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Mozambique, data.civsoc.Mozambique), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Mozambique, data.civsoc.Mozambique), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Mozambique, synth.res = synth.civsoc.Mozambique, tr.intake = 1975, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('NICARAGUA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Nicaragua, data.civsoc.Nicaragua), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Nicaragua, data.civsoc.Nicaragua), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Nicaragua, synth.res = synth.civsoc.Nicaragua, tr.intake = 1979, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('RWANDA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Rwanda, data.civsoc.Rwanda), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Rwanda, data.civsoc.Rwanda), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Rwanda, synth.res = synth.civsoc.Rwanda, tr.intake = 1994, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('U.S.S.R.', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Russia, data.civsoc.Russia), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Russia, data.civsoc.Russia), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Russia, synth.res = synth.civsoc.Russia, tr.intake = 1917, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('VIETNAM', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Vietnam, data.civsoc.Vietnam), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Vietnam, data.civsoc.Vietnam), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Vietnam, synth.res = synth.civsoc.Vietnam, tr.intake = 1954, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
  
  title_af = paste('YUGOSLAVIA 1945', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.civsoc.Yugoslavia, data.civsoc.Yugoslavia), 3), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.civsoc.Yugoslavia, data.civsoc.Yugoslavia), 3), sep = '')
  path.plot.cus(dataprep.res = data.civsoc.Yugoslavia, synth.res = synth.civsoc.Yugoslavia, tr.intake = 1945, Ylim = c(0, 1), Main = title_af, Ylab = 'Core Civil Society Index')
  abline(v=0, lty=2)
dev.off()

